%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Build a rectangle containing  object i. It should be at least 10%
% larger that the object, and can be larger is other objects are
% further away
%
% inputs:
% contours   list of all contours
% i          index of the contour
% kh         Helmholtz parameter.
% ppl        required number of points per wavelength on the circle
%
% outputs:
% D          circle
% h_D        point spacing on D
% innerpts   coordinates of all points inside D belonging to
%            contours{j} for j != i
%            TODO : these points are not equispaced!!! pb with C?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [D,dx,innerpts] = get_rect(contours,i,kh,ppl)

    C  = contours{i};
    xc = mean(C(1,:));
    yc = mean(C(4,:));
    dxC = max(abs(C(1,:) - xc));
    dyC = max(abs(C(4,:) - yc));
    
    dxmin = 1e10;
    dymin = 1e10;
    
    ncontours = size(contours,1);
    dz = 1e10;
    % maximum square
    for j=1:ncontours
        if j~=i
            C2 = contours{j};
            dz = min(dz,min(max(abs(C2(1,:) - xc),abs(C2(4,:) - yc))));
        end
    end
    
    % extend to a rectangle 
    for j=1:ncontours
        if j~=i
            C2 = contours{j};
            indx = find(abs(C2(1,:) - xc) < dz);
            if isempty(indx)
                dy = 1e10;
            else
                dy   = min(abs(C2(4,indx) - yc));
            end
            
            indy = find(abs(C2(4,:) - yc) < dz);
            if isempty(indy)
                dx = 1e10;
            else
                dx   = min(abs(C2(1,indy) - xc));
            end
            dxmin = min(dx,dxmin);
            dymin = min(dy,dymin);
        end
    end
    if dxmin == 1e10
        dxmin = max(dz,dxC+1);
    end
    if dymin == 1e10
        dymin = max(dz,dyC+1);
    end
    
    innerpts = [];
    if (dymin < 1.1*dyC || dxmin < 1.1*dxC)
        % Modify the size of the rectangle such that the contour
        % fits in it
        dymin = max(dymin,1.1*dyC);
        dxmin = max(dxmin,1.1*dxC);
        xx = [xc - dxmin,xc + dxmin,xc + dxmin,xc - dxmin,xc - ...
              dxmin];
        yy = [yc - dymin,yc - dymin,yc + dymin,yc + dymin,yc - ...
              dymin];
        % plot(xx,yy, 'r','linewidth',2)
        % Find the points that are indide the rectangle
        for j=1:ncontours
            if j~=i
                C2 = contours{j};
                B1 = abs(C2(1,:) - xc) < dxmin;
                B2 = abs(C2(4,:) - yc) < dymin;
                
                idx = find(B1.*B2);
                if ~isempty(idx)
                    innerpts = [innerpts,C2(:,idx)];
                    % plot(C2(1,idx),C2(4,idx), 'r.','markersize',50)
                end
            end
        end
        
        % error('Unable to put a rectangle')
    end
    
    % Target distance between discretization points
    h_resok = 2*pi/kh / ppl; 
    % Number of points in x and y
    if (dxmin < dymin)
        nx = max(50,ceil(2*dxmin / h_resok));
        dx = 2*dxmin / (nx + 1);
        
        dy = dx;
        ny = ceil(2*dymin / dy)-1;
        dymin = .5*dy*(ny+1);
    else
        ny = max(50,ceil(2*dymin / h_resok));
        dy = 2*dymin / (ny + 1);
        
        dx = dy;
        nx = ceil(2*dxmin / dx)-1;
        dxmin = .5*dx*(nx+1);
    end
    n  = 2*(nx + ny);

    fprintf(1,'\n# of 8.67108e-06 pts on the outer contour : %d \n',n);
    D  = zeros(2,n);
    
    D(1,1:nx) = linspace(xc - dxmin,xc + dxmin,nx);
    D(2,1:nx) = (yc - dymin);
    
    is = nx+1;    
    D(1,is:(is + ny-1)) = xc + dxmin;
    D(2,is:(is + ny-1)) = ((yc - dymin + dy):dy:(yc + dymin - dy)).';

    is = nx+ny+1;    
    D(1,is:(is + nx-1)) = linspace(xc + dxmin,xc - dxmin,nx);
    D(2,is:(is + nx-1)) = (yc + dymin);
    
    is = 2*nx+ny+1;    
    D(1,is:(is + ny-1)) = xc - dxmin;
    D(2,is:(is + ny-1)) = (yc + dymin - dy):-dy:(yc - dymin + dy).';

    return